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Abstract. We investigate the dynamical instabilities of a resonator coupled to a 
superconducting single-electron transistor (SSET) tuned to the Josephson quasiparticle 
(JQP) resonance. Starting from the quantum master equation of the system, we use a 
standard semiclassical approximation to derive a closed set of mean field equations 
which describe the average dynamics of the resonator and SSET charge. Using 
amplitude and phase coordinates for the resonator and assuming that the amplitude 
changes much more slowly than the phase, we explore the instabilities which arise in 
the resonator dynamics as a function of coupling to the SSET, detuning from the JQP 
resonance and the resonator frequency. We find that the locations (in parameter space) 
and sizes of the limit cycle states predicted by the mean field equations agree well with 
numerical solutions of the full master equation for sufficiently weak SSET-resonator 
coupling. The mean field equations also give a good qualitative description of the set of 
dynamical transitions in the resonator state that occur as the coupling is progressively 
increased. 
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1. Introduction 

Nanoelectromechanical systems in which a high frequency mechanical resonator is 
coupled to a mesoscopic conductor p] have been predicted to display a wide variety 
of different dynamical behaviours depending on the nature of the conductor. When a 
mechanical resonator is linearly coupled to the transport electrons in either a quantum 
point contact [21 [3j H] , or a normal state single-electron transistor [3 El [7], the electrons 
act on the resonator like an effective thermal bath [2, El El El El H] . Under such 
circumstances, the resonator is damped and reaches an almost Gaussian steady state 
whose width is set by the fluctuations in the motion of the charges through the 
conductor. In contrast, where the electro-mechanical coupling is non- linear [91 [10 j, or 
the conductor is close to a transport resonance [T2l [13], the mechanical resonator 
can be driven by the electrons into states of self-sustaining oscillation. 

In this paper we analyze the instabilities that arise in the dynamics of a mechanical 
resonator coupled to a superconducting single-electron transistor [HI [HI EE! EE] (SSET) 
operated in the vicinity of a transport resonance. In the SSET, transport resonances 
occur when states of the SSET island differing by one Cooper pair are degenerate so 
that coherent Cooper pair tunnelling between the island and one of the leads is possible. 
The simplest such resonance, which we concentrate on here, is called the Josephson 
quasiparticle ( JQP) resonance pH [181 El H3] an d involves a cycle of processes in which 
current flows via a combination of coherent Josephson tunnelling between the SSET 
island and one of the leads, followed by incoherent quasiparticle decays into the other 
junction. 

When a SSET tuned to the vicinity of the JQP resonance is coupled to a mechanical 
resonator, the resonator dynamics is very sensitive to the precise choice of bias point 
for the SSET [HUH]. F° r operating points on one side of the JQP resonance, the SSET 
damps the resonator and can be regarded as an effective thermal bath, behaviour which 
was confirmed in recent experiments [16] . In contrast, for operating points on the other 
side of the resonance the electrical degrees of freedom pump the resonator leading to the 
possibility of states of self-sustaining oscillation. In this regime the resonator dynamics 
can be investigated either by numerical solution of the quantum master equation [12] , or 
provided the resonator is much slower than the electrical degrees of freedom (as was the 
case in recent experiments [HIT6]), an effective Fokker-Planck equation for the resonator 
can be derived [HJ [13] . 

Numerical solution of the SSET-resonator master equation [12] revealed interesting 
similarities with a quantum optical device known as the micromaser. In a 
micromaser [T9l [20] , a cavity resonator is driven by the passage of a steady stream 
of excited two-level atoms. As the atom-cavity interaction is increased, the resonator 
undergoes a sequence of dynamical transitions leading in some cases to non-classical 
steady states. A corresponding set of dynamical transitions was found to occur in the 
SSET-resonator system together with regions of non-classicality. 

Here we use an alternative approach, namely a mean field description to analyze 
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the instabilities in the dynamics of the SSET-resonator system. This kind of approach 
has been used extensively in the analysis of non-linear quantum optical systems where 
its usefulness has been well established [2lJ 22] • The mean field equations we derive 
provide a relatively compact description of the system and their stability properties are 
readily analyzed using the techniques of classical dynamical systems theory. Although 
the mean field description neglects some of the correlations in the system, comparison 
with numerical results from the full master equation show that for sufficiently weak 
coupling the mean field theory is close to being quantitatively correct. Although 
quantitative agreement is poor at stronger couplings, the mean field equations still give 
a good qualitative description of the system's dynamics displaying a similar sequence of 
transitions to that found numerically [12J. 

Although we will use language appropriate to a nanomechanical resonator 
throughout this paper, we anticipate that much of our analysis will also be relevant to 
the case of a superconducting resonator coupled to a SSET. Photon assisted satellites 
of JQP peaks have been observed in experiments in which a SSET was coupled to 
microstrip transmission line [23]. Furthermore, recent experiments have demonstrated 
coherent coupling between a superconducting stripline resonator and a Cooper-pair box 
with a coupling Hamiltonian between the resonator and the charges on the box very 
similar to the mechanical case [23] . Such systems would only need to be modified to 
allow quasiparticle tunnelling off the Cooper-pair box [25] to become analogous to the 
device considered here and hence it seems likely that a range of dynamical instabilities 
similar to that which we find for the mechanical case could also be produced in a 
superconducting resonator. 

The outline of this paper is as follows. In Sec. 2 we introduce our model of the 
SSET-resonator system and give the appropriate quantum mechanical master equation. 
The mean field equations are derived in Sec. 3. In Sec. 4 we transform the mean 
field equations into plane polar coordinates and exploit the fact that the amplitude of 
the resonator oscillations is slowly changing to derive an effective amplitude dependent 
damping of the resonator due to the SSET. We then show that this quantity can be used 
to predict the presence of limit cycles in the resonator dynamics. We compare the limit 
cycle solutions predicted by the mean field theory with numerical calculations using the 
full master equation in Sec. 5. Then in Sec. 6 we briefly consider the implications of 
our theoretical calculations for experiments on nanomechanical- SSET systems. Finally, 
in Sec. 7 we draw our conclusions. In the Appendixes we give additional details about 
the derivation of the master equation for the SSET-resonator system and the stability 
analysis of the mean field equations. 

2. Master Equation 

The SSET-resonator system we consider is shown schematically in figure UK. The 
SSET consists of an island linked to leads by tunnel junctions with resistances Rj and 
capacitances Cj which are taken to be equal for simplicity. The mechanical resonator is 
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Figure 1. (a) Circuit diagram for the coupled SET-resonator system. Under the 
bias voltage shown the electrons flow through the SSET island from left to right, (b) 
Summary of charge processes involved in the JQP cycle. 



treated as a single-mode harmonic oscillator, with frequency u, which forms a movable 
gate capacitor. The equilibrium position of the resonator is a distance d from the 
SSET island and we assume that the displacement of the resonator with respect to 
this equilibrium position, x, is always small in comparison (i.e. \x\ <C d). Under these 
circumstances the gate capacitance can be expanded up to just linear order and we have 
C g {x) — C g (l — x / d) , where C g is the capacitance at x = 0. 

The JQP resonance [T71 [T8| [TT| [T4"] which we are interested in here occurs when 
two conditions are met. Firstly, at the centre of the resonance there is no change in the 
energy of the system when a Cooper pair tunnels between one of the leads and the SSET 
island. Secondly, the bias voltage, Vd s , must be large enough to allow quasiparticle decay 
processes between the island and the other lead to occur. The charge processes involved 
in the JQP resonance we consider are summarized in figure [lb. Josephson tunnelling 
across the left junction leads to coherent oscillations between the SSET island states |0) 
and 1 2) which differ by a single Cooper pair. These oscillations are interrupted by the 
tunnelling of a quasiparticle from the island into the right lead which takes the island 
into charge state |1). A further quasiparticle then tunnels from the island to the right 
lead, returning the island to state |0) and the cycle begins again. The large electrostatic 
charging energy of the small SSET island and the carefully controlled bias voltage ensure 
that at low temperatures all other charge processes are strongly suppressed. 

The master equation for the SSET island charge and resonator is obtained from the 
full Hamiltonian of the system by tracing out the quasiparticle degrees of freedom |17j . 
The steps involved in this derivation (together with details of the approximations and 
simplifications involved) are sketched in Appendix A. The final result is a master 
equation of the form [12] 



The evolution consists of a coherent part, described by the Hamiltonian, H ccn together 
with two dissipative terms £i ea ds and ^damping which describe quasiparticle decay from 
the island and the surroundings of the resonator respectively. The effective Hamiltonian 




(1) 
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is given by [T2] 

tf co = A£|2>(2|-^(|0><2| + |2>(0|) (2) 
v 2 1 

+ + rmwV + mcu 2 x s x (|1)(1| + 2|2)(2|) , 
2m 2 

where AE is the energy difference between states |2) and |0), Ej is the Josephson energy 
of the superconductor, x and p are the canonical position and momentum operators of 
the resonator and m is the effective resonator mass. The resonator-SSET coupling is 
described by the length-scale x s which measures the shift in the equilibrium position of 
the resonator brought about by adding a single electronic charge to the SSET island [Hj. 
The tunnelling of quasiparticles from the island is described by 

Clea d sP= -^[{|2)(2| + |l)(l|,p} + (3) 

-2(|l)(2| + |0)(l|)p(|2)(l| + |l)(0|)], 

where V is the quasiparticle decay rate. Note that this is a simplified expression which 
uses a single rate for the two decay processes and neglects both the position dependence 
of the quasiparticle rates and their dependence on bias point (these approximations are 
discussed in more detail in Appendix A). Simplified in this way, the master equation 
is also essentially equivalent to that used to describe a resonator coupled to a double 
quantum dot [26] . 

The term which describes the effect of the resonator's surroundings is 

CdampingP = _ "7^f[ X ' />}+] ~ Te ^ (w + l/2)[g, [x, p]], (4) 

where ^ ex t represents the external damping and n the equilibrium average resonator 
occupation that would occur in the absence of the SSET. We use this form of the 
oscillator damping kernel as opposed to the quantum optical form which we have used 
elsewhere [12] as although it is less convenient for numerical calculations, it leads to the 
correct (translationally invariant) classical limit, which is essential if we wish to derive 
appropriate mean field equations [271 12H] • 

The behaviour of the resonator depends very sensitively on the sign of AE. 
For AE > 0, Cooper pairs tend to absorb energy from the resonator damping its 
motion [29j [TI], [14] . In contrast, for AE < and when the resonator is slow (i.e. w < T) 
the resonator tends to absorb energy from the Cooper pairs and it is in this regime 
that instabilities can occur. For faster resonator speeds [12J (uj/T > 1), absorbtion of 
energy by the resonator from the SSET, and hence the location of instabilities, becomes 
concentrated around points where AE ~ —jhuo with j an integer. 

3. Mean Field Equations 

The mean field equations for the SSET-resonator system consist of the set of equations 
of motion for the expectation values of all the relevant SSET and resonator operators. 
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They are derived by multiplying the master equation by each operator in turn and 
taking the trace: 



where the SSET charge operators are defined as pij = \i)(j\ with i,j = 0, 1,2 and the 
expectation values are defined by (...)= Tr[p . . .]. Note that, as discussed in Appendix 
A, the pi2,Pio components of the density operator decouple from the rest and can safely 
be neglected from the mean field equations (and also the master equation itself) as they 
do not affect the resonator. 

In contrast to the simpler case of a resonator coupled to a normal state SET [6, 30J, 
the equations of motion of the first moments do not form a closed set. In other words, 
the dynamics of the first moments depend on the behaviour of higher moments leading 
to an infinite hierarchy of equations of motion for progressively higher order moments 
of the system operators. In order to derive a simple set of dynamical equations we need 
to make a semiclassical approximation whereby the expectation value of a product of 
two operators is replaced by the corresponding product of the expectation values of the 
individual operators [211 122] . Thus in this case we substitute (x)(po 2 ) for (xp 02 ). This 
approximation cannot be justified rigorously. Indeed, dropping the correlations between 
x and P02 means that we lose the ability to describe the noise in the system (which is 
determined by the behaviour of the higher moments). Nevertheless, the semiclassical 
approximation is well-known as a way of deriving a set of dynamical equations which 
typically capture many of the important elements of the dynamics of the corresponding 
quantum system [22]. In this case, we find that the much simpler mean field equations 
which result from the semiclassical approximation provide a very useful qualitative and, 
for sufficiently weak SSET-resonator coupling, quantitative description of the different 
dynamical transitions which the resonator can undergo. 

Thus, having made the semiclassical approximation, and using the conservation of 
probability ((poo) + (pn) + (P22) = 1) to eliminate one of the charge variables, we obtain 
the following closed set of mean field equations 




v 



x 



V 



J 1 (X + X s [pu + 2p 22 }) - JextV 



(11) 

(12) 




(13) 
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P =±(AE + hw^a-^(3 + ^(2p 22 + Pll -l) (14) 

Pll = r(p22 - Pll) (15) 
E 

P22 = ~ - Tp 2 2, (16) 



where x 2 q = H/(2mu) and we have dropped the angled brackets for convenience. The 
quantities a and (3 are defined as the real and imaginary parts of (p 02 ) respectively. 

Despite the decoupling of the second moment in the above equations, the mean 
field equations reproduce many of the features found in the dynamics of the full master 
equation. The time evolution of these equations reveals the resonator relaxing to a 
fixed point for some values of the parameters. For other values, the oscillations grow at 
first, before settling into a limit cycle, i.e. an oscillation at fixed amplitude. The mean 
field equations also show clear evidence of bistability for certain parameter values: in 
these cases the long-time behaviour of the resonator depends on the choice of initial 
conditions. 

The fixed point solution of the mean field equations is obtained by setting the time 
derivatives to zero in equations (|TT|) - (fl6|) . As we discuss in Appendix B, the stability of 
the fixed point can be established using standard techniques [22J. 

4. Analysis in radial coordinates 

Although a straightforward stability analysis of the mean field equations can tell us 
quite a lot about the dynamics of the system, we find that transforming the mean field 
equations into plane-polar coordinates and making one further simplifying assumption 
allows us to proceed much further. The assumption we make is that the SSET- 
resonator coupling and external damping are sufficiently weak that the resonator's 
energy changes much more slowly than either its phase or the charge state of the 
SSET, conditions which are readily met for most practical implementations of the SSET- 
resonator system. This type of approximation has already proved useful in describing a 
variety of nanoelectromechanical and optomechanical systems |10[ [TJ1 ED] . In terms of 
our analysis of the mean field equations, the assumption of a wide separation of time 
scales between the evolution of the amplitude and phase of the resonator allows us to 
derive an effective equation of motion for the resonator amplitude from which we can 
determine the number and location of stable limit cycle solutions [31] . 

We proceed by rewriting the mean field equations [equations (TXT]) - (TT6]) ] in plane 
polar co-ordinates (A, 0) which describe the amplitude and phase of the resonator, 
defined through the relations x — Xf p = Acos<p and v = uAsmcj), where Xf p is the fixed 
point resonator position. In terms of (A, <p) the mean field equations take the form, 

A = - 7 exf Asin 2 - u[x fp + x s (p n + 2p 22 )} sin0 (17) 
<j> = - u - 7 ex . t sin cos - -^[x fp + x s (p u + 2p 22 )} cos0 (18) 

PXl = r (p 2 2 - Pll) (19) 
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(20) 
(21) 

(22) 



Note that Xf p is equal to — x s (pn + 2p 22 ) when pu and p 22 have their fixed point values. 

We now assume that the evolution of the resonator amplitude, A, is sufficiently 
slow that during a single period of oscillation it can be treated as a constant and 
that the phase evolution over the same time can be approximated by <fi = out. The 
dynamics of the electronic degrees of freedom can then be obtained for this constant 
amplitude and steadily evolving phase. The resulting forced dynamics of the SSET 
charge variables are then averaged over the resonator period to calculate their effect on 
the resonator amplitude A which we characterize by an amplitude dependent damping 
term p2JJ [111 [31] , Jsset{A). This slow- A approximation will certainly be appropriate in 
the vicinity of a limit cycle solution of the mean field equations and will be valid more 
generally so long as the free (uncoupled) evolution of both the resonator and charge 
degrees of freedom is much faster than the rate of change of the resonator amplitude in 
the coupled system, i.e. for u>,T 3> ^ e xt and sufficiently weak SSET-resonator coupling. 
Whilst there is no simple way of evaluating the strength of the SSET resonator coupling 
at which this approximation breaks down, from equations ( fTTj) and ({TBI , we expect 
that this approach will be valid for a given amplitude, A, if x s <C A (note that if the 
populations pu and p 22 remain much less than unity then the conditions on x s will 
be less restrictive). Furthermore, it is clear that in order to be consistent with these 
assumptions we should obtain an effective damping •yssEriA) whose magnitude is much 
less than u and Y. 



4-1. Solving the electronic dynamics 

With a fixed amplitude, A, and a harmonically oscillating phase, = out, the electronic 
degrees of freedom form a set of four coupled differential equations with time dependent 
coefficients. The dynamics of even this simplified system is still non-trivial, but we can 
make progress by making use of our assumption that the electronic degrees of freedom 
relax rapidly compared to the resonator amplitude. We assume that transients in the 
charge dynamics can be neglected and hence that the effect on the amplitude of the 
resonator is dominated by the periodic response of the charge degrees of freedom to the 
harmonic drive. 

To extract the relevant periodic solutions for the charges we rewrite the equations 
of motion in terms of a Fourier series consisting of harmonics of the resonator frequency, 
defined by, e.g. 

+00 

Pn(t) = £ P n ne iwt - 
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The resulting equations for the Fourier coefficients of the electronic variables are, 

iamjft = r {p n 22 - p n n ) (23) 
iunp n 22 = - Yp n 22 - ^f3 n (24) 

iumoT = - + F - ^(P n+1 + /T" 1 ) - \a n (25) 

*onF = + + §(2^ 2 + Pn - 5 nfi ) 

+ u^(a^ + a^)-^. (26) 

By solving for p™i,p 22 , ct n in terms of (3 n , we can rewrite these equations so that we have 
an equation for f3 n in terms of (3 n+1 , /3 n_1 , /3 n+2 , f3 n ~ 2 . This is equivalent to a matrix 
equation involving a band-diagonal matrix with five non-zero diagonals. The matrix 
equation we must solve is, 

= M(3_- Ejd/(2h), (27) 

where (3 is the vector of the coefficients (3 n (with n running from — oo to +oo), M 
represents the matrix of terms derived from equations (I2"3"|) - (j2"6"|) . and d is a vector 
representing the Kronecker delta, 5 H; o, i.e. d(0) = 1 and d(n ^ 0) = 0. 

Solving the Fourier series is equivalent to inverting the matrix M and the coefficients 
(3 n are given by, 

P = E J M- 1 d/(2h). (28) 

The components decay rapidly towards zero for \n\ 3> 1 and so the matrix can be 
truncated and solved numerically with little error. In practice we found that calculating 
the coefficients up to n = ±80 proved more than adequate for the parameters we 
considered. 

Examples of the oscillations in the SSET charge driven by the resonator are shown 
in figure [2] for the case where the resonator frequency matches the quasiparticle decay 
rate. We use a single simplified set of SSET parameters in all numerical calculations 
to illustrate the behaviour of our model [H] EG]: T = Vd s /eRj, Rj = h/e 2 and 
Ej = hVd s /(16eRj). The SSET-resonator coupling strength is parameterized by the 
dimensionless quantity k = muj 2 x 2 /eVd s - 

The charge oscillations in figure [2] resemble the behaviour of a periodically kicked 

damped oscillator, but more specifically, at larger amplitudes they are very similar 

to an oscillator whose frequency is time-dependent. Analogous oscillations have been 

seen in the amplitude of radiation within an optomechanical cavity [32| |3T] . We can 

establish the unforced behaviour of the charge oscillations for a given amplitude using 

equations (fl9l) - (l2"2"]) and treating the phase as constant (0 = 0). The unforced SSET 

charge oscillations have a frequency which increases linearly with the amplitude, A. For 

sufficiently small Ej, we can approximate the frequency by 

AE ojx s , . „ . , 

u>e*-jr- + 1 Tr{xf P + A) (29) 
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Figure 2. Time dependent oscillations of the SSET charge described by 2p22 + Pn, 
calculated numerically using the Fourier series solution (solid curve) and using the 
analytic approximation described in Sec. 14.31 (dashed curve). A harmonic oscillation 
with the resonator frequency is plotted (as a dotted curve) for comparison. The 
parameters used were lu/T = l,AE/eV ds = -0.1, K = 0.04, with A/x s =2.4 (i), 9.4 
(ii), 21.8 (hi) and 40.7 (iv). 



and the decay rate is ~ T/2. Thus as the amplitude of the resonator increases, the 
number of oscillations in the SSET charge degrees of freedom during each resonator 
period increases accordingly. 



4-2. Effective damping 

Having calculated the response of the SSET charge to the periodic driving provided by 
the resonator, we can now calculate how the resonator, in turn, responds to the dynamics 
of the SSET. Our assumption is that the change in the amplitude of the resonator 
oscillations is much slower than the oscillations themselves. We can therefore average 
the effect of the SSET charge dynamics on the resonator over a single resonator period, 
over which time the amplitude of the resonator motion can be treated as constant. 
Integrating equation ( TT71) over the resonator period gives, 



dA _ ~f ext A 
dt ~ 2 



^- J dt'cux s (pn + 2p 2 2)smujt', 

(30) 



Dynamical instabilities of a resonator 

0.10 
0.08 
^"0.06 
K 0.04 
0.02 
0.00 



11 




10 20^^30 40 50 

s 



Figure 3. Behaviour of the SSET effective damping, 'jsset, as a function of the 
resonator amplitude, A/x s . The curves show —"/ssetA/x s (in units where T = 1) 
calculated numerically (full curve) and an analytic approximation (dashed curve) 
calculated using equation The parameters used are the same as in figure [5] 

Also plotted are lines indicating the external damping for jext/F = 0.005 and 0.001 
with the latter being the less steep. Intersections between the curves and one of the 
lines indicates the presence of a limit cycle. The points labeled (i)-(iv) correspond to 
the oscillations shown in figure El 



where the bar on the amplitude, A, indicates that the equation is only valid on time 
scales longer than 2n/uj. In terms of the Fourier series, the integration eliminates all of 
the terms except those with n = ±1, hence the equation of motion for A can be written 

as 

^ = - ?f A - x^Im [ P \ X {A) + 2pl 2 (A)} . (31) 

The second term can be interpreted as an amplitude dependent damping arising from 
the interaction with the SSET [I0| [TT], ^ S set(A)A = 2uJx s hn[p\ x {A) + 2p\ 2 {A)] as 
it appears on the same footing as the external damping rate, ^ext- However, when 
1sset{A) < it means that the charges transfer energy to the resonator. The resonator 
has a constant amplitude (but not phase) whenever dA/dt = 0, and hence supports 
limit cycle solutions whenever this condition is met for A ^ 0. 

Figure [3] shows [33] jsset(A)A calculated numerically as a function of A, along 
with various values of j e xt- Limit cycle solutions exist for the values of A where the 
curve of —Jsset(A)A is crossed by a given line representing 7 ea; tA The stability of the 
limit cycle solutions depends on the gradient of dA/dt with respect to A in the usual 
way, hence a limit cycle is stable where this gradient is negative [34] . Thus we can see by 
inspection from figure [3] that for the case where 7 ex t/r = 0.001, there are three possible 
limit cycle solutions two of which are stable (the largest and smallest amplitude ones). 

It is clear from figure [3] that it is the oscillations in , Jsset{A) as a function of A 
that lead to the existence of more than one limit cycle solution [35] for sufficiently weak 
j ex t. The oscillations in 'Jsset(A-) can in turn be understood as arising from changes in 
the commensurability of the SSET charge oscillations with the resonator frequency as 
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the amplitude is increased. This effect can be seen clearly by comparing the oscillations 
in figure [2] with those in figure [3j At the first peak in —JssetA, the charge undergoes 
one oscillation during the resonator period, a number which increases by two at each 
subsequent peak in — 'JssetA. 



4-3. Approximate analytic solution 

Although the SSET charge is bounded within a narrow range of values, its oscillations (as 
shown in figure [2]) nevertheless resemble those of a driven harmonic system. Essentially 
this is because the parameters chosen are such that the Josephson energy is not sufficient 
to allow the charge to saturate over the time scale of the oscillations. This similarity 
suggests that it should be possible to find an approximate solution by reducing equations 
(fi~9i)-(f2"2"l) to an equation for an appropriately driven harmonic oscillator. We can then 
follow the approach used in [31] where the harmonic dynamics of an optical cavity 
coupled to an oscillating mirror was analyzed. 

For an uncoupled SSET [TT| [14"] . it is straightforward to show that in the limit 
where Ej/hT <C 1 the populations of the upper charge states, pn and P22 always remain 
much less than unity. Furthermore, examining the numerical evolution of the mean 
field equations we find that this remains true in the coupled system. This suggests that 
we can simplify the analysis in the limit where Ej/hT <C 1 by neglecting the pu and 
P22 terms in equation fl22l . Making this approximation, and again assuming that the 
resonator damping is much slower than the other time scales, we find that equations 
( T2~T]) and (1221) reduce to the equations of motion for a damped harmonic oscillator with 
a time dependent frequency term. Using P02 = ct + i/3, we obtain a single complex 
equation of motion, 

n .Ej 



P02 



AE uox s . . 
+ —2~\ x fp + Acosuit) 



h x 2 



2 J yo 2h 



(32) 

which can be solved by use of a Fourier series, P02 — £ %e & lu}nt Po2 where the tilde 
indicates that the Fourier series includes a global phase shift 9{t) = z sin uot and 
z = Ax s /x 2 q . The value of Pq 2 can be written in terms of Bessel functions of the first 
kind, pQ 2 = ^ n Jn{,—z\ where the parameter ip n is given by, 

^ = iEj/n (33) 

2{iun — i(AE/h + u)x s Xf p /x 2 q ) — T/2) 
As before, the damping is calculated from the Fourier coefficients of pn and P22, 
which in turn can be calculated from f} n using equations (1231) and ([21 



Ajsset = -uJX s —flm 
n 



1 V 



(34) 



T + iu (T + iu) 2 ^ 

In order to make use of this Bessel function solution, we must convert between Fourier 
series with and without global phase shifts. We find, 

^ = Ji E {^ mj m+n(z) - *J m - n {z)) Jm{z). (35) 
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This approximate solution can be used to calculate Pn{t),P22(t) (figure [2]) and hence the 
effective damping (figure [3]). It is noticeable in figure [3] that even for our relatively large 
choice of Ej/hT ~ 0.4, the analytical approximation agrees well with the numerics at 
large amplitudes (i.e. for A/x s ^> 1). This is because as the amplitude of the resonator 
increases, the driven oscillations which develop in the SSET charge degrees of freedom 
become progressively faster and the populations pn(t) and become ever smaller 

as can be seen in figure [2j Using smaller values of Ej leads rapidly to better agreement 
at small amplitudes. 

5. Comparison with numerical results 

Much of the usefulness of our analysis of the mean field equations rests on the degree to 
which this simplified description of the SSET-resonator system faithfully reproduces the 
behaviour seen in the full numerical solution of the master equation |12j . The mean field 
equations allow us to calculate the amplitude of limit cycles in the resonator dynamics 
as a function of all the various parameters of the system, to determine which of them are 
stable, and calculate the associated SSET current. Of course the mean field description 
does not describe the noise in the SSET-resonator system and hence it cannot tell us 
the degree to which a particular limit cycle may be occupied in regions of parameter 
space where there is more than one stable limit cycle (or a stable fixed point solution 
coexists with one stable limit cycle). 

5.1. Size of limit cycles 

We begin by considering sufficiently weak couplings (for a given value of 7 e:rf ) that the 
resonator is limited to at most a single limit cycle state [12] and examine when such 
states develop and their sizes as a function of the detuning from resonance, AE, and 
the relative speed of the resonator, u/T. In figure H] the size of the stable limit cycles 
calculated using the mean field equations [i.e., using equation (131]) ] is compared with the 
numerical solution of the full master equation as a function of AE for various resonator 
frequencies. 

The numerical solution of the master equation gives us the full steady-state density 
matrix of the system, p ss , in terms of which the probability of finding the resonator in a 
given Fock state, \n), is simply P{n) = Tr[\n)(n\p ss ]. Because of the ensemble average 
and the presence of phase noise, limit cycles do not appear as periodic features in the 
dynamics of the density matrix, but they can be identified as peaks in the distribution 
P{n) above the ground state energy (i.e. n > 0). [These peaks typically correspond 
to individual rings in the associated Wigner function representation of the density 
matrix [12].] We have used two quantities from the P(n) distribution to compare with 
the mean field results: the average number of resonator quanta, (n) = ^2 n nP(n), and 
the n-value of the peak in the distribution, which we define as ni c . In order to compare 
these values with the stable limit cycle amplitudes calculated in the mean field theory, 
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Figure 4. The size of the limit cycles as a function of AE/eVd s for a range of values 
of u/T, with, 7ext/r = 0.002, k = 0.01. The prediction of the mean field theory, 
n rn f is compared with ni c and (n) obtained from the numerical solution of the master 
equation. The darker (blue) shaded regions indicate limit cycle states and the lighter 
(yellow) shaded region indicates the regions of bistability, as determined from the 
master equation. Also shown for comparison are the regions where a stable fixed 
point solution exists calculated using the fixed point analysis described in Appendix 
B (arrows). 



we express the latter in terms of resonator quanta, n m f = A 2 {muj/2K). 

It is clear from figure H] that the mean field equations prove to be a rather good 
predictor of the locations and sizes of the limit cycles in the full system dynamics as 
can be seen by the relatively good agreement between n m f and ni c . It is interesting to 
note that a comparison of the sizes of the limit cycles with (n) works well in regions 
where the resonator energy distribution P(n) is relatively concentrated. Thus n m f is 
quite close to (n) when there is a single, large, limit cycle state present, but not in the 
vicinity of the transitions where the limit cycles begin to form or in regions where a 
limit cycle state coexists with a stable fixed point (bistable regions). 

Despite the generally good agreement between n mp and ni c , figure H] reveals that 
there are small differences in the locations of the values of AE at which the limit cycles 
appear (or disappear) predicted by the mean field equations and the master equation. 
This can partly explained by the fact that our analysis is most appropriate when the 
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Figure 5. The size of the limit cycles as a function of k 1 / 2 calculated numerically and 
using the mean field equations, with u>/T = 1, 7 er /r = 0.001, AE/eVd s = —0.1. The 
solid green line denotes the mean field solution, and the dotted blue curve indicates the 
size of the limit cycles calculated using the approximate Bessel function solution. The 
background shading shows the numerically calculated steady state distribution of the 
resonator, P{n) 1 with dark red indicating highest probability. The white dots indicate 
the peaks in P{n). 



limit cycles are large on a scale set by x s (as discussed in Sec. H]). However, the differences 
in onset points for the limit cycles exist not just in cases where the limit cycles grow 
continuously from zero (around AE = 0), but also when they emerge with a relatively 
large radius. An interesting possible explanation for these observations is that the exact 
location of the dynamical transitions may in fact depend quite sensitively on the noise 
in the system — an element which is of course missing in our mean field analysis [22l [13] . 

We now turn to compare the mean field analysis with numerical predictions for 
a given detuning, AE, as the SSET-resonator coupling is increased. For resonator 
frequencies of the order of the quasiparticle decay rate, numerical solution of the master 
equation showed that increasing the coupling could lead to a sequence of transitions 
marked by the appearance of increasing numbers of peaks in the steady-state distribution 

P{n)m- 

In figure [5] we plot the size of the stable limit cycles calculated both numerically 
using the Fourier series solution of the mean field equations and using the Bessel 
function expression [equation (1341) ] as a function of k 1 ^ 2 for uj/Y = 1. The predictions 
of the mean field equations are compared with the locations of the peaks in the 
numerically calculated resonator distribution, P(n). The mean field calculation shows 
good qualitative agreement with the full numerics, showing a series of bifurcations and 
multiple limit cycles as the coupling is increased. For weaker coupling, k 1//2 < 0.1, the 
mean field calculation accurately predicts the size of the limit cycle. Notice, however, 
that the Bessel function approximation for the limit cycle differs appreciably from the 
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full mean field for the first limit cycle, but otherwise matches the mean field result 
closely. As before, reducing the value of Ej improves the accuracy of the Bessel function 
approximation. 

The appearance of successive stable limit cycle solutions as the coupling is increased 
in figure [5] is readily understood in terms of the analysis of the charge oscillations given 
in Sec. 14.11 and the associated oscillations in 'Jsset(A) (illustrated in figure [3]). From 
equation (1291) we see that increasing the SSET-resonator coupling (i.e. increasing x s ) 
increases the frequency of the charge oscillations thus changing their commensurability 
with the mechanical period. This effectively compresses the oscillations of '')sset(A) as 
a function of A (i.e. they occur with a progressively smaller period measured in terms 
of amplitude). Thus for fixed jext, increasing the coupling means that more and more 
stable limit cycle solutions occur and those already present move to smaller and smaller 
sizes. 

5.2. SSET Current 

It is also possible to calculate the current through the SSET using the mean field 
equations. The current is generated by quasiparticle tunnelling out of the states |1) 
and |2), which leads to a time dependent tunnel current across the right junction, 
I(t) = T(pn(t) +P22{t))- Therefore, when the resonator is in a limit cycle state of a 
particular amplitude, the corresponding oscillations inp n (t) andp 2 2(^) [see figure[2] will 
be passed on to the tunnel current. The average current (defined as either an average 
over one period of mechanical oscillation, or over an ensemble of systems) of course 
will not reflect these oscillations, but it will depend on the amplitude of the resonator's 
motion. For a given resonator amplitude, the corresponding current is given by the 
Fourier coefficients, (I)/e = T^p^A) +£>2 2 (^))- 

The average current calculated using the mean field equations only agrees well 
with that calculated using the full master equation when the steady state of the latter 
predicts a narrow width to the distribution P(n). However, we expect that the current 
calculated within the mean field picture will be useful in understanding the current noise 
in the SSET [36]. For example, we might expect signatures of the different frequencies 
of charge oscillations corresponding to different resonator states to appear in the current 
noise spectrum. 

6. Practical Implications 

The mean field analysis can be used to provide simple estimates of the kinds of dynamical 
instabilities that may be seen in a particular experiment. Recent experiments on a 
22MHz nanomechanical resonator coupled to a SSET were carried out [16] in the regime 
where u/T <C 1. Previous calculations [TH [13], [H] showed that it should be possible 
to see the transition from the fixed point to a limit cycle state in such systems and 
a similar conclusion is reached using the mean field equations. Indeed some evidence 
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Figure 6. Behaviour of the SSET effective damping, jsset, as a function of the 
resonator amplitude, A/x s . The parameters used are as described in the text below. 
The curve shows —jssetA/x s , while the straight lines indicate the external damping 
for 7 eK t/r = 2.25 x 10~ 5 and 5 x 1CT 6 with the latter being the less steep. 



for instabilities in the resonator motion was found in these experiments, although the 
primary focus of the work was on the stable regimes [TBI [37] . 

It is natural to use the mean field approach to estimate whether the regime 
where the resonator dynamics involves more than one limit cycle solution is accessible 
experimentally For the nanomechanical resonator and SSET parameters in the 
experiments of Naik et al. [16] we find that the wide separation of time-scales, u/Y -C 1, 
means that the values of the SSET-resonator coupling and 7^ necessary to reach this 
regime are well beyond currently achievable values. However, the interaction between the 
SSET and resonator is much stronger when uj/T ~ 1, which suggests that the observation 
of more than one limit cycle may be possible when the resonator and quasiparticle decay 
timescales are more closely matched. 

Although nanomechanical resonators with frequencies ~ 1GHz have been 
produced, [38] integration with SET electronics has only been achieved for resonators 
up to ~ 100MHz [7j. However, it is possible to reduce the quasiparticle tunnelling rate, 
T, substantially by increasing the resistance of the relevant tunnel junction and rates 
< 2 x 10 8 s _1 have been demonstrated [25] . Therefore, as an example, we consider a 
nanomechanical resonator with fundamental frequency u/2ir = 100MHz coupled to a 
SSET in which all the charge processes are rather slower than usual. Quasiparticle 
tunnelling in the SSET is assumed to occur across a junction with very high resistance, 
5MQ, and a Josephson energy at the other junction which is tuned (e.g. using 
the method employed by Nakamura et al. [25]) to be Ej = 1 x 10~ 3 eVd. s . For the 
other parameters of the system we choose values which are within the general range 
explored in recent experiments [TBI EES]. The other SSET parameters we assume are 
Ec = 175yueV, eVd s = 700/ieV, and A = 200yueV. For the resonator we assume a mass 
m = 6.8 x 10~ 16 kg, with a SSET-resonator separation d = lOOnm, a SSET-resonator 



Dynamical instabilities of a resonator 



18 



capacitance C g = 34aF and coupling voltage V g = IV. We note that this choice of 
parameters takes us outside the limit Ej/hT <C 1 and we must rely on the numerical 
solution to the full Fourier series (rather than use the Bessel function approximation). 

In figure [6] we plot the effective SSET damping as a function of resonator amplitude 
for a detuning AE = — 1 x 10~ 3 eVd. s . From figure 6, we see that we first get two stable 
limit cycles around r y e xt/T — 2.5 x 10~ 4 , corresponding to a resonator quality factor of 
2.9 x 10 4 , which should be accessible experimentally. The observation of multiple limit 
cycles of course requires that the fluctuations in resonator amplitude are not so large as 
to wash out the difference in amplitudes between the cycles. Although a full calculation 
of the noise in the system is beyond the mean field theory as presented, we can make an 
estimate of the length scale of the fluctuations due to external thermal noise which, for 
the parameters in figure [6] and a temperature of 30mK we find to be about 8 A ~ 136x s 
i.e. about 10% of the separation in amplitude of the limit cycles. This suggests that 
multistability will occur within a region accessible by experiment though the conditions 
required to see it are quite demanding. 

7. Conclusions 

The mean field analysis presented here provides a simplified description of the SSET- 
resonator system. In particular, the mean field approach provides a good description of 
the dynamical instabilities which the resonator can undergo. For relatively weak SSET- 
resonator couplings the mean field equations describe the onset and size of the first limit 
cycle state quite accurately for a wide range of resonator frequencies. As the coupling 
is increased, the mean field equations predict the emergence of a succession of limit 
cycle states of different sizes. Numerical calculations based on the full master equation 
reveal that although the mean field equations become progressively less accurate as 
the coupling is increased they still give a good qualitative description of the dynamical 
behaviour. Furthermore, the mean field analysis allows us to relate the appearance 
of further limit cycles to changes in the commensurability of oscillations in the SSET 
charge with the period of the resonator. 
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Appendix A. SSET-resonator master equation 

In this appendix we outline some of the steps in the derivation of the master equation, 
equation (CQ) and the approximations made. Essentially the same master equation was 
introduced in Ref. [14J and used again in [12], but details of its derivation have not yet 
been given. 
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In essence the calculation is a straightforward generalization of the treatment of 
the JQP cycle given in Ref. |T7] to include coupling to a resonator. The Hamiltonian 
for the SSET-resonator system can be written as the sum of several parts, 

H = H c + H 3 + H R + H (lp + H T (A.l) 

which correspond to the charging energy of the SSET, He, coherent Josephson coupling 
between the left lead and the island, Hj, the energy of the resonator, Hr, the energy 
of the quasiparticles in the leads and island of the SSET, H qp and the quasiparticle 
tunnelling Hamiltonian which couples the leads and the island, H^. 

As discussed in Sec. [2j for relatively small resonator displacements the resonator- 
charge coupling can be linearized and hence the charging energy of the SSET written 
as [10] 

H c = i Ec ( N2 ~~ 2NN 9) - neVds + mu 2 x s xN] \N)(N\ <g> \n)(n\ (A.2) 

N,n 

where E c = e 2 /(2Cj + C g ), N g = (C g V g + CjV ds )/e and x s = 2E c C g V g /(emuj 2 d). The 
macroscopic charge variables N and n, which correspond to charge states \N) and \n), 
are the number of excess electronic charges on the SSET island and the number of 
electrons which have tunnelled off the island via the right hand junction respectively. 
The Hamiltonian of the resonator is simply that of a harmonic oscillator, 

Hr = f- + \rmo 2 x 2 . (A.3) 
2m 2 

For the bias configuration we consider, Josephson tunnelling is only important between 
the left lead and the island as tunnelling between the right lead and the island involves 
energy differences ~ 2eVd s 3> Ej because of the bias voltage applied (we neglect the 
possibility of the resonator having large enough energy to assist in overcoming this 
barrier). Thus we include just the coherent tunnelling between the island the left lead 
(i.e. coherent tunnelling changes N but not n), 

H 3 = -J2^(\N)(N + 2\ + \N + 2)(N\) . (A.4) 

TV 

The quasiparticle Hamiltonian is given by [T7] 

-^qp = tkac\ aa Ckaoi (A-5) 

a=L,R,I k,<7 

where the sum a runs over the three pieces of superconductor (left lead L, right lead R, 
and island I) and the operator c\ aL creates a quasiparticle in the left lead with energy 
ekL, momentum k and spin state a. Quasiparticle tunnelling between the island and the 
leads is described by the Hamiltonian, 

£ (e'^^Xj + e+^xty (A.6) 



H^ = 

j=L,R 



where 



Xj = T k q c\ ja c qIa (A.7) 

k,q,a 
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creates a quasiparticle in state k in lead j and destroys one in state q in the island with 
tunnelling amplitude T^ q . The operators e ±l ^ j ^ 2 describe the associated change in the 
macroscopic charge variables, in terms of which they are written 

6 t#h/2 = ^ \ n ± l)(n| ® |JV =f l)(iV| (A.8) 



AT,n 



= J2\Nt1)(N\. (A.9) 

N 

The master equation for the macroscopic island charge N and resonator is 
obtained by tracing out the quasiparticle degrees of freedom and the count variable 
n. The procedure is very similar to that used to derive the analogous master 
equation for a normal-state SET [30]- The derivation uses the standard Born-Markov 
approach [21] which involves the assumptions that the tunnelling Hamiltonian Ht is a 
weak perturbation on the quasiparticle Hamiltonian and that the quasiparticles relax 
back to their unperturbed distributions after tunnelling faster than any other time-scale 
in the problem. 

Further simplification is achieved by limiting the analysis to include only those 
states which are accessible to the system in the zero temperature limit. The charging 
energy of the SSET island, E c is typically much larger than Ej, hence the Josephson 
tunnelling can be neglected for all states except the two charge states which are selected 
at the JQP resonance (by tuning [JT] N g , see equation IA. 21) to be almost degenerate. For 
simplicity we consider states |0) and |2), which differ by one Cooper pair, corresponding 
to resonance (for a fixed gate) at N g — 1. For quasiparticle decay to occur, the energy 
gained when a particle tunnels from the island to a lead must be > 2 A, where A is 
the superconducting gap [14] . At the JQP resonance the voltage Vd s is chosen so that 
only two processes are allowed: tunnelling from the island into the right lead between 
states |2) and |1), and between |1) and |0). The displacement of the resonator produces 
changes in the electrostatic energy differences involved in quasiparticle tunnelling which 
are assumed to be small with respect to the values for a fixed gate. This is the reason 
for our choice of a dimensionless of coupling, k = muo 2 x 2 s / eVd s <C 1, which measures the 
energy associated with coupling to the resonator in terms of the typical energy scale 
associated with the unperturbed quasiparticle tunnelling rates. The question of when 
other charge states become accessible (i.e. for very strong SSET-resonator couplings or 
high enough temperatures), and what effect they have on the dynamics is an interesting 
one, but we do not consider it here. 

Taking all these factors into account, we obtain a master equation confined to the 
space of the three charge states |2), |1) and |0), 

p=~[H m ,p] (A.10) 



T{E 2t 



2 



X ± [{|2)(2|,p} + - |l)(2|p(|l)(0| + |2)(1|) - (|0)(1| + |l)(2|)p|2)(l|] 



m^r'(£ 2 ,i) [ {x | 2 )(2|,p} + - x|l)(2|p(|l)(0| + |2)(1|) - (|0)(1| + |l)(2|)p|2)(l|x] 
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2 

where T(E) is the quasiparticle tunnel rate for energy E (this is the energy gained by 
the quasiparticle when it tunnels from island to lead) with the relevant energies for 
quasiparticle tunnelling from |2) to |1) and from |1) to |0) given by |14j. 

E 2 ,i = eV ds + 2£ c (3/2 - N g ) (A. 11) 

£i,o = eV ds + 2E c (l/2 - N g ), (A.12) 

respectively. We have assumed that the position dependent changes in the energy 
differences are small with respect to the typical energy scale eVd s , so that we can expand 
the tunnelling rates to first order [T3] , 

T(E + muj 2 x s x) ~ T(E) + mu 2 x s xT'(E) (A. 13) 

where T'(E) = dT(E)/dE. The Hamiltonian in equation flA.101) takes the simplified 
form 

H co = AE\2)(2\ + £ x |l)(l| - ^ (|2)(0| + |0)(2|) 

+ mu?x s x [|1)(1| + 2|2)(2|] + H R , (A.14) 



where 



AE = 4E C (1 — N g ) (A. 15) 

Ex =AE c (l-2N g ). (A.16) 

The master equation we have obtained [equation flA.lOj) ] can be analyzed 
numerically quite easily. However, it is useful to make further simplifications which make 
it easier to identify the essential physics of the system. The off-diagonal elements of the 
density matrix (0|p|l), (l|p|2) (together with their complex conjugates) decouple from 
the dynamics of the other parts of the density matrix and can be dropped, together with 
the Ex term in H co , as they play no role in the resonator dynamics. We also assume that 
the difference between the two quasiparticle decay rates and their variations with bias 
point can be neglected. Whilst the difference between these rates (and their dependence 
on the bias point) may be important in the analysis of a given experimental system, it 
is not essential to our theoretical analysis which seeks to describe the basic features of 
the system in the simplest possible way. 

Finally, we also drop the position dependent parts of the quasiparticle tunnelling 
terms. All the important dynamics arises from the position dependence in the coherent 
part of the master equation which modulates the detuning from resonance AE (which 
can be arbitrarily small). In contrast, the position dependent coupling which appears 
in the quasiparticle decay terms is expected in practice to be much weaker than the 
coherent position coupling and provide only a small modulation of the quasiparticle 
rates [TTJ [HJ 13] . Although numerically (within the master equation formalism) we find 
that including the position dependence of the quasiparticle rates can eventually lead 
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to quite large changes in the steady state distribution function for sufficiently large 
k, the general pattern and range of resonator behaviours (including the existence of 
regions where the resonator state is number-squeezed) remains essentially the same. 
Furthermore, corrections arising from the position dependence of the quasiparticle rates 
become progressively smaller as Rj/(h/e 2 ) is increased [TT1 fT3] . Note, however, that 
the parameters used in the main text (and in [12]) are chosen to best illustrate the 
behaviour of our simplified model and we have not attempted to in addition minimize 
the corrections that would arise if the position dependence of the quasiparticle rates 
were included. 



Appendix B. Stability analysis of the fixed point 

In this appendix we describe the analysis of the fixed point solutions of the mean field 
equations [equations (|TTT)-(|T6l)]. The fixed point solution, j3f p , is found from the solution 
of the following cubic equation, 



wj \vj T Hfp t - w w yfp 

( AE 2 3E 2 T r\ „ Ej 

in terms of which the fixed point values of the other dynamical variables are readily 
obtained. In particular, the fixed point resonator position is given by, 
3E 

x fp = -^-x s f3 fp . (B.2) 

It is obvious from equation ( IB. II) that there is only a single solution for (3 in the limit 
n — > 0. We have checked numerically that there remains only one physically acceptable 
solution over the range of parameters studied here. 

In order to establish the stability of the fixed point solution we need to calculate 
the Jacobian matrix of the mean field equations at the fixed point |34j . The stability of 
the system is then determined by the eigenvalues of this matrix, A, which are solutions 
of the characteristic equation which can be written as 

A 6 + a 5 A 5 + a 4 A 4 + a 3 A 3 + a 2 A 2 + aiA 1 + a = (B.3) 

The coefficients do, a±, etc. are functions of the various parameters of the system 
determined by taking the determinant of the appropriate Jacobian. Whilst it is possible 
to determine the stability of the system by simply calculating all the eigenvalues 
numerically it turns out that we can establish the stability of the fixed point solution 
using a somewhat simpler approach. 

Assuming we start from a stable region and slowly change one of the parameters 
of the system, a limit cycle of the resonator begins to form when a pair of complex 
eigenvalues cross the real axis (Hopf bifurcation) [3U [22J, H2]. The locations of these 
transitions can be determined by using the fact that at these points the characteristic 
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equation must support a solution of the form A = iv (with v a positive real number). 
This requirement can be re-expressed in terms of the condition on the corefficients |22[|4"2] 

/ a\ + a 5 (a a 3 - aia 2 ) \ 2 / a\ + a 5 (a a 3 - a x a 2 ) \ _ 

/ = «5 —7 r— - o 3 I —7 — + ai = (B.4) 

yas^aoas — a\a±) + aia% J \(i5(,ao a 5 — 0104) + 0103 / 

The regions of stability marked by arrows in figure H] were determined by evaluating / for 
each set of parameters. In stable regions / < 0, changing sign at the Hopf bifurcations. 
Notice that our analysis is centred on the fixed point solution; it tells us nothing about 
the possible coexistence of stable fixed point and limit cycle solutions and hence does 
not describe the regions of bistability in figure HI 
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